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We theoretically study the development of quantum turbulence from two counter-propagating su- 
perfluids of miscible Bose-Einstein condensates by numerically solving the coupled Gross-Pitaevskii 
equations. When the relative velocity exceeds a critical value, the countersuperflow becomes unsta- 
ble and quantized vortices are nucleated, which leads to isotropic quantum turbulence consisting of 
two superflows. It is shown that the binary turbulence can be realized experimentally in a trapped 
system. 
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Hydrodynamic instability develops into turbulence or 
turbulent flow due to the complicated dynamics of eddies 
of various length scales [1]. This phenomenon appears in 
diverse fields of physics, such as magneto- hydrodynamic 
fluids in plasma physics [2], relativistic fluids in as- 
trophysics J3| , and quantum fluids in low temperature 
physics 0, [J Turbulence has been studied thoroughly in 
classical fluid dynamics as the most important unsolved 
problem of classical physics. Recently, quantum turbu- 
lence (QT), turbulence in supernuids, has attracted inter- 
est in low temperature physics [6]. QT consists of topo- 
logical defects, vortices with definite quantized circula- 
tion, and can provide a simple prototype for understand- 
ing classical turbulence. QT has been historically studied 
in superfluid 4 He and 3 He and has been numerically stud- 
ied in atomic Bose-Einstein condensates (BECs) @, Q 
and has recently been realized in experiments 

Early studies on QT were devoted primarily to ther- 
mal counterflow in superfluid 4 He, where the hydrody- 
namics is usually described using a two-fluid model in 
which the system consists of an inviscid superfluid com- 
ponent and a viscous normal fluid component [To| . Ther- 
mal counterflow is driven by an injected heat current, 
where the two components flow in opposite directions. 
For a sufficiently large relative velocity between the two 
components, the laminar counterflow state becomes un- 
stable, developing into a superfluid turbulent state. In 
the transition, remnant vortices, quantized vortices at- 
tached to the container wall, are stretched by the mutual 
friction between the superfluid and normal fluid compo- 
nents, growing into a tangle through reconnect ions with 
other vortices 11,0]. 

We consider a similar case in two-component BECs, 
QT arising from the instability of the counterflow 
of two superfluids, namely, countersuperflow instabil- 
ity. Counter-superflow instability has been theoreti- 
cally studied in several multicomponent miscible super- 
fluid systems, such as helium superfluids [l3l4l6j, mixture 
BECs of cold ato ms Il7j , and nucleon superflows in rotat- 
ing neutron stars [18|. Atomic BEC systems have great 
advantages over other superfluid systems, because mod- 



ern experimental techniques enable us to easily control 
superflows and directly visualize the dynamics of topo- 
logical defects therein. Recently, Hamner et al. experi- 
mentally realized countersuperflow instability for the first 
time, and observed that shocks and dark-bright solitons 
nucleate via the instability in quasi-one-dimensional two- 
component BECs [l9]. We expect that countersuperflow 
instability exhibits novel hydrodynamic instability as was 
reported in @, [H[ , and causes quantum turbulence in two- 
component BECs as in the case of the thermal counter- 
flow turbulence. In the context of the research of study 
QT, it is remarkable that the superfluid velocity field in 
two-component BECs is continuous in contrast to that in 
single-component superfluids [io| . The continuous veloc- 
ity field enable us to calculate the important statistical 
quantities in classical turbulence such as enstrophy and 
helicity, which are not well defined in single-component 
superfluid systems. Hence we expect that countersuper- 
flow instability in two-component BECs will open new 
directions for research into unique turbulent states in 
multicomponent superfluids. 

In this work we investigate the development of coun- 
tersuperflow instability into QT and its decay in two- 
component BECs at zero temperature by numerically 
solving the coupled Gross-Pitaevskii (GP) equations. In 
contrast to QT in 4 He thermal counter flow, the insta- 
bility develops into QT without remnant vortices in our 
system. In the countersuperflow instability, quantized 
vortices are intrinsically nucleated in the bulk and then 
stretched by the momentum exchange between the two 
condensates [Fig. [T] (a-d)]. Eventually, quantized vor- 
tices of the two condensates get tangled with each other 
forming a binary QT [Fig. [U^e)]. By calculating the time 
evolution of enstrophy of the binary QT, we investigate 
its isotropy and the similarity with classical turbulence. 
We also demonstrate that the binary QT arising from 
the countersuperflow instability can be realized by the 
current experimental technique in cold atoms. 

Two-component BECs at low temperatures are de- 
scribed in mean- field theory [2l| by the macroscopic wave 
functions Vj(t,r) = y/rij^r)^ 6 ^^ (j = 1,2) with the 
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FIG. 1: (Color online) Nonlinear dynamics of countersuperflow instability. The figures show the isosurfaces of the density 
ni = 0.05no of the 1st component. The instability causes (b) disk-shaped density pulses which (c) transform into vortex rings. 
The vortex rings are enlarged by the momentum exchange between the two condensates, and then (d) the vortex interactions 
and reconnections deform the rings, leading to (e) binary quantum turbulence, (f) Then the vortex tangle decays. The 
isosurfaces with n<i — 0.05no of the 2nd component is also plotted in (e) in a different color. We solved the GP equations with 
the relative velocity U12 = 4.71 x ^Jgno/m on 64 3 grid. The box size is (11 x ^) 3 with the characteristic vortex core size 
fo = ft I ^/mno(g + #12). Time t is normalized as t' = gnot/h. 



particle densities rij and phases 0j, where the index j 
refers to the jth component. The wave functions are 
governed by the coupled GP equations 
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Here rrij and Vj(r) are the mass and the external po- 
tential of the jth component. The coefficient gjk = 
27rh 2 ajk/rrijk represents the atomic interaction with 



and the s-wave scattering length a 



jk 



between the jth and kth components. Our analysis satis- 
fies the conditions #11(722 > fl?9 and gjj > that two mis- 
cible condensates are stable [2l|. Note that, since there is 
no exchange of particles between the two components, the 
equation of continuity of mass is satisfied for each com- 
0, where we use the mass density 

Pj v j 



V Jo 



rrijUj and the mass flux Jj 
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jth component. 

Countersuperflow instability. Before discussing 
trapped systems, we first elucidate the countersuperflow 
instability in an isolated uniform system with Vj = 0. 
Because of Galilean invariance, we may neglect the trans- 
lational motion of the whole system without loss of gener- 
ality: P\ + P2 =0 with Pj = J Jjdr. The wave function 
of the jth component in a stationary state may then be 
n o e ji(mjUj-r-njt) w jth the chemical 
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potential fij and the superfluid velocity U\ 

The countersuperflow state is dynamically unstable 
when the magnitude of the relative velocity U12 = \U\ — 
Z7 2 | exceeds a critical value U c . The critical relative ve- 
locity of the countersuperflow instability is essentially dif- 
ferent from the critical velocity of the Landau instability 
[lOj, at which the frictionless flow becomes thermody- 
namically unstable. In the Landau instability, a super- 
fluid exchanges momentum with a "rigid" environment 
(e.g., container wall or normal component), and the su- 
perfluid system is treated in the framework of a canonical 
ensemble. In contrast, the countersuperflow instability is 
purely an "internal" instability of the isolated system 
of two superfluids without any influence from the exter- 



nal environment. Thus the countersuperflow instability 
should lead to an exchange of momentum between the 
two components to reduce their relative motion, keeping 
the total momentum conserved in the uniform system. 

The momentum exchange must arise from a gradual 
excitation of internal motions, that is, from the appear- 
ance of excitations in the superfluids. The momentum 
exchange due to the excitations is investigated by solving 
the Bogoliubov-de Gennes (BdG) equations, which are 
obtained by linearizing the GP Eqs. (pQ) with respect to a 
collective excitation SVj = e^ m ^' r -^ {u j e i< < q " r - UJ ^ - 
v * e -i(q-r-ujt)^ ar0 und the stationary state. The total 
momentum density S J carried by the excitation can be 
defined as S J — 8J\ + 8J2 with the change of the mo- 
mentum density SJj = hq(\uj\ 2 — \vj\ 2 ) of the jth com- 
ponent. The unstable modes, which trigger the counter- 
superflow instability, should satisfy the condition S J = 
with SJi = — 5J2 7^ due to the law of momentum 
conservation. Therefore, since the the countersuperflow 
instability is the dynamic instability triggered by the un- 
stable modes with the imaginary part Im uj ^ [13], 
the amplification of the unstable modes exponentially en- 
hance the momentum exchange. 

Nonlinear dynamics of the transition towards binary 
QT. We demonstrate the characteristic nonlinear dy- 
namics of the countersuperflow instability in an uniform 
system by numerically solving the GP Eqs. (pQ) with 
Vj = in a periodic box. The parameters are set to 
be m = mi = ra 2 , g = gu #22, and g 12 = 0.9g, 
similar to the experimental situation discussed later. In 
the numerical simulation, the initial state is prepared 
by adding a small random noise to the wave function 
y°j(t 0) with n = n? = n° 2 and U x = -U 2 . Figure 
[U shows the time development of low-density isosurfaces 
with 77,1 (r) = 0.05no, where U± = ^Ui2X with a unit 
vector x along the x-axis. Because of the symmetric pa- 
rameter setting, the 2nd component behaves similarly to 
the 1st component (not shown in Fig. [T]). 

In the early stage of the dynamics, the amplification 
of the unstable modes creates disk-shaped low-density re- 
gions that orientate in the x direction [Fig. [T](b)]. The 
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lowest density inside the disk region reaches zero, cre- 
ating a local dark soliton [l9| . The soliton in the jth 
component transforms into a vortex ring [22[ with a mo- 
mentum antiparallel to the initial velocity Uj [Fig. H(c)]. 
The distribution of the vortex rings can be determined 
from the dispersion relation uj(q) obtained by diagonal- 
izing the BdG equations; 



hW = e 2 + -tfqp^ ± y/e*h*q*U* 2 + 4£ 2 2 n 2 e 2 , (2) 

where we use e 2 = e (e + 2#n ), e = |^<? 2 , and 
q 2 = qf\ + Q± with qf|| = \q • x\ and q± > 0. The instabil- 
ity condition a; 2 < is reduced to e_ < hq\\Ui2/2 < e + 
with e± = yeo[eo + 2(g ± #12)^0)] • This inequation de- 
termines the unstable region in the wave number space 
(q\\iQ±) when C/12 exceeds the critical relative velocity 
U c = U- with U± = 2^/1 ± gi2/g\/9no/m. The unsta- 
ble modes have characteristic wave numbers ^|| < 7^12 
and ^ i^^ 12 ^ or l ar S e relative velocity U12 > U+. 
The length scale ~ with ft = 27rh/m then character- 
izes both the radii of the vortex rings and the intervals 
between the rings along U12 just after the vortex nucle- 
ations. Thus, the vortex line density l v after the insta- 
bility is roughly estimated to be l v ~ £/ 2 2 /ft 2 , which is 
controllable by changing the relative velocity U\2- 

The momentum exchange accelerates after the vortex 
ring nucleation [Fig. 12(a)]. The motion of the vortices 
then dominates the exchange. Since the momentum car- 
ried by a vortex ring increases with its radius, the radii 
of the nucleated vortex rings increase with time for the 
momentum exchange. This dynamics is similar to that of 
quantized vortices under the thermal counterflow of he- 
lium where the vortices are dragged by the mutual 
friction between the superfluid and normal fluid compo- 
nents. Our numerical simulations show that the "mu- 
tual friction" between the two condensates drags the vor- 
tices of each condensate under the counter superflow. The 
"mutual friction" may be understood as an interaction 
between a vortex line and excitations. 

When the vortex rings become large, the interaction 
between the vortex rings deforms the rings and vortex 
reconnections occur, which depresses the momentum ex- 
change [FigJT] (d) and Fig 12(a)]. These effects make the 
vortex dynamics much complicated, leading to binary QT 
where vortices of both components are tangled with each 
other [Fig. H(e)]. Then the momentum exchange almost 
stops and each component has about zero momentum on 
average. 

Statistical properties of the binary QT. The interac- 
tion between vortices of different components, which is 
repulsive for gi2 > [20], becomes important when the 
mean intervortex distance is small. In fully developed 
binary QT, where the vortices of the two components 
are strongly correlated, it is plausible to regard the two 
condensates as a single fluid. We then introduce the to- 
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FIG. 2: (a) Time evolution of momentum (Pj) x — 
y/gmno(Pj) x along the x-axis of the jth component. The 
two black arrows show the time around which the vortex 
nucleations (VN) and reconnections (VR) start, (b) Nu- 
merical plots of velocity \v\ — \Jgno/m\v'\ and densities 
p — mnop\ pi = mnop'i, pi — mn^p'2 (inset) as a functions 
of the distance r = hr / \fgmno from the core of the quasi- 
Rankine vortex in a stationary state, (c) Decay of anisotropic 
parameters Q\\ = qoQu, Q± = qoQ± and enstrophy Q = qoQ f 
with qo — (gno/h) 2 (inset). 



tal mass density p = mini + rri2n2 and the total super- 
fluid mass-current velocity v = (Ji + J2)/ P- Since the 
low-density regions along vortex cores of one component 
are filled with the other components due to the repul- 
sive interaction between the components, the total mass 
density p is leveled. Consequently, the velocity field v 
around the coreless vortex mimics that of the Rankine 
vortex where the velocity field v is not irrotational 
and thus the vorticity field uj = V x v is continuous [Fig. 
lib)]. 

It is convenient to investigate the statistical properties 
of the binary QT using the quantity Q = ^7 f oj 2 dr with 
the system volume V, which is called the enstrophy in 
classical turbulence Analogizing with the discussion 
of helium QT in thermal counterflow , we define the 
anisotropic parameters 
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where Q = Qu + Q±. If the tangle is isotropic, we have 
Q\\ = Q_l/2. Conversely if the tangle consists of curves 
lying only in planes perpendicular to Z7 12 , then Qu = 
and Q± = Q. The time evolution of Q\\ and Qj_ in 
the countersuperflow QT are plotted in Fig. [2(c). The 
value Q± first appears with the vortex ring nucleation 
and increases due to the vortex line extension by the 
mutual friction keeping Q\\ = 0. The vortex line mixing 
due to the vortex interactions and reconnections increases 
Q\\ to ~ Q_l/2 when the turbulence becomes isotropic. 
Then the turbulence starts to decay with the enstrophy 
Q keeping its isotropy Q\\ ~ Q_\_/2. 

We can directly analogize the binary QT with classi- 
cal turbulence by investigating the time evolution of the 
decaying enstrophy. In a model for the decay of classical 



turbulence in a finite system [23], enstrophy decays like 
Q(t) oc t -11 / 5 in the early stage of the decay process. 
The enstrophy Q of the binary QT follows this decay in 
the early stage [inset of Fig. EJc)], which means that the 
binary QT consisting of the quasi-Rankine vortices can 
mimic the classical turbulence as a turbulence of a single 
fluid. Such a viewpoint has not formerly been proposed 
and presents another direction for the analogy between 
classical turbulence and QT. To confirm this viewpoint, 
we have to find further similarities in different statistical 
uantities i.e. energy spectrum [7[, velocity distribution 
, etc. A detailed analysis of binary QT will be reported 
elsewhere. 

Counter- super flow QT in a trapped system. Counter- 
superflow QT can be realized experimentally by employ- 
ing a technique similar to that of Ref. [19| , where relative 
motion of the two condensates is induced by the Zee- 
man shift. We consider two-component BECs of 87 Rb in 
the hyperfine states |F, mj?) = and |2, — 1), which 

can satisfy the miscible condition gnQ22 > <Zi2 by us- 



ing the intercomponent Feshbach resonance [24[ . We use 



an isotropic harmonic trap Vj(t = 0,r) = yc^r 2 with 
uj r = 2tt x 8Hz and set the particle number to N = 10 6 
for both components and the scattering lengths to a = 
an = CL22 = 5.2 nm and a\2 — 0.9a. Since the full three- 
dimensional numerical simulation takes too much time 
and computer memory, we calculated a two-dimensional 
simulation considering the z = cross section of the 
condensates to investigate qualitatively whether the QT 
transition is observable experimentally. 

The condensates are first located at the minimum of 
the trap [Fig. [3ja)]. A magnetic field gradient —151 
mG/cm is applied at t = ms, which accelerates the two 
condensates in opposite directions along the x-axis due 
to Zeeman shifts. Although countersuperflow instability 
can occur in the region where the two condensates over- 
lap, this collision process is not sufficient for the insta- 
bility to grow into QT [Fig. [3jb)]. However, this process 
is important for introducing sufficient seeds that develop 
into QT through the next collision process, described be- 
low. The field gradient is turned off at t = 11.7 ms 
and the two condensates again return to the potential 
center r = and collide there. Then, the countersu- 
perflow instability forms a quasiperiodic density pattern 
consisting of many dark solitons [Fig. 02(c)]. The soli- 
tons transform into vortex pairs, which corresponds to 
the vortex ring nucleations in three-dimensional space. 
The nonlinear vortex interaction leads to a turbulent 
state [Fig. 03(d)] even if there is initially very small noise 
Srij(r) ~ 10 _4 nj(r). Our simulation shows that the 
whole dynamic of the transition to binary QT from coun- 
tersuperflow instability can be observed experimentally. 
A detailed analysis of trapped systems will be reported 
elsewhere. 
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Scientific Research from JSPS (Grant No. 21340104). 




FIG. 3: (Color online) Time evolution of the density profile 
n[ of the 1st component in the transition to binary QT from 
countersuperflow instability in a trapped system. The black 
lines represent the density contour n' 2 = 0.004 of the 2nd com- 
ponent. The density nj is normalized as nj — n~ a) 
with R T = (15Na/a h ) 1/5 . The field of view is 29.2a,, x 20.8a^ 
with ah — \/h/ (mujr). 
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